%%MATLAB Code for:
%%Successive Over-Relaxation (SOR) 
%%Solution of Laplace's Equation
%%on Rectangular Region
%======================================
%%Instructor: Dr. Joe LoVetri
%%Developer: Mohammad Shafieipour
%======================================

%Dirichlet conditions all around rectangular region
%Calculation of total flux emanating from region

%clear MATLAB environment
clc;
clear all;
close all;

%declareations
phi=ones(11,11);%solution matrix
epsilon=10^-16;%relative displacement norm
omega=1.5;%acceleration factor (1<omega<2)
max_iter=10^2;%maximum iterations

%inpose Dirichlet conditions
for i=1:11
    %keep in mind that in MATLAB indexing starts from 1
    phi(i,1)=0;%bottom side
    phi(i,11)=100;%top side
    phi(1,i)=0;%left side
    phi(11,i)=0;%right side
end

%initial guess
%set all interior potential values to 1
%since phi is created using ones() function, it is already done

%count number of iterations
iter=1;
Err_norm=99999;

while iter<max_iter && Err_norm>epsilon
    phi_norm=0;%solution norm
    delta_norm=0;%displacement norm
    for i=2:10
        for j=2:10
            Residual=(phi(i-1,j)+phi(i+1,j)+phi(i,j+1)+phi(i,j-1))/4;
            Acc_residual=omega*(Residual-phi(i,j));
            phi(i,j)=phi(i,j)+Acc_residual;
            phi_norm=phi_norm+abs(phi(i,j));%calculate current solution norm
            delta_norm=delta_norm+abs(Acc_residual);%calculate displacement norm 
        end
    end
    Err_norm=delta_norm./phi_norm;%calculate relative displacement norm
    iter=iter+1;
end
contourf(phi.',20)
colormap (flipud(colormap(gray)))
colorbar
%calculate total flux emanating from region
flux=0;
for i=2:10
    flux=flux+phi(1,i)+phi(i,1)+phi(11,i)+phi(i,11);%outside potential
    flux=flux-phi(2,i)-phi(i,2)-phi(10,i)-phi(i,10);%inside potential
end

%display the output
disp ('The calculated solution is:');
disp (phi);
fprintf('Total flux emanating from region = %d\n',flux);







